Yukawa bosons in two-dimensional harmonic confinement 
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The ground state property of Yukawa Bose fluid confined in a radial harmonic trap is studied. The 
calculation were carried out using the Density Functional Theory formalism within the Kohn-Sham 
scheme. The excess correlation energy for this inhomegeneous fluid is approximated via the Local 
Density Approximation. A comparison is also made with the Gross-Piteavskii model. We found 
that the system of bosons interacting in term of Yukawa potential in a harmonic trap is energetically 
favourable compared to the ones interacting via contact delta potential. 

PACS numbers: 61.20.Ja, 67.40.Db, 31.15.Ew 



YUKAWA BOSE FLUID, AN 
INTRODUCTION 



A system of N Bose particles interacting via the 
Yukawa potential is called Yukawa Bose fluid (for short 
YBF). For three dimensional-space D = 3 the Yukawa 
potential takes the form Vs(f)=(ea/f) exp(— r/a), where 
e is an energy scale and a is a length scale having the 
meaning of a screening length. Whereas for the two- 
dimensional case however, it is described by V^r) = 
eK (r/a) where K (x) is the modified Bcssel function 
that decays as exp(— x) / y/i at large distance while it be- 
haves in a repulsive manner as — ln(x) at short distances. 
In other word the Yukawa potential presents a combina- 
tion of short range with soft-core type potential. 

Early numerical calculation (quanta! Monte Carlo sim- 
ulations) in studying the ground state property of the 
three-dimensional YKF can be traced back to the work 
of Ceperley and coworkers [H, 0]. There were a spurt 
of interest in numerical calculation of 2D ground state 
property of the 2D- YBF in the 90's following the idea of 
Nelson and Seung Q who have shown that the statistical 
mechanics of the flux-line lattice (FLL) of high-T c super- 
conductors can be studied through an appropriate map- 
ping onto the 2D- YBF. Magro and Ceperley 0] have per- 
formed a Diffusion Monte Carlo (DMC) and Variational 
Monte Carlo (VMC) numerical simulation to calibrate 
the ground state vortex properties and phases (liquid or 
solid) that correspond to the density and kinetic energy 
of the system. Similar work related to the first order 
phase transition of Abrikosov lattice to liquid of vortices 
have been reported by Nordborg and Blatter [j| using 
the path integral Monte Carlo method. In the late 90's 
other methods such as the STLS model @, 0] appeared 
in studying the ground state properties of the 2D- YBF. 
It has to be noted that all these works were spiralling 
around the system of homogeneous fluid. 

Our work will focus on the harmonically trapped sys- 
tem extending the work based on the quantal Monte 
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FIG. 1: The phase diagram of Yukawa bosons reproduced 
from Magro and Ceperley Q- The pluses are transition points 
computed with DMC. The dashed line at high densities is the 
scaling law: A ~ 0.04/y/p. 



Carlo calculations of Magro and Ceperley [1] . Of partic- 
ular interest to us is the phase diagram obtained by them 
for the parameters (p, A) at transition points displayed 
in Fig. JT]). Here A indicates the Dc Boer dimension- 
less parameter defined by A 2 = ti 2 /2ma 2 e while p and m 
denote the reduced density and atomic mass respectively. 

Magro and Ceperley [J] observed that for a particular 
region A > 0.09 the system is dominated by kinetic en- 
ergy and does not crystallize. Below this threshold how- 
ever a peculiar behaviour of reentrant liquid has been ob- 
served. Which means that system is in the liquid phase 
at very low density as well at high density. The crys- 
tal melts on compression and expansion. It is known 
that at low density the system is in liquid phase due to 
the fact that the screened potential cannot bind to solid 
phase particles that are on average too far apart. For 
the high density the crystal is known to melt similar to 
the mechanism of thermal melting of a Wigner crystal. 
However in this work we are dealing with the liquid phase 
of Yukawa bosons only. Our calculation rely on the data 
of the liquid phase near the transition point from the 
phase diagram Fig. ([T]) as an input to the Local Density 
Approximation (LDA). The use of Density Functional 
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theory (DFT)-LDA in this work is only a matter of com- 
putational convenience primarily due to the availability 
of the exact Monte-Carlo [H data of the homogeneous 
system. Infact to our knowledge this is the first attempt 
to study the YBF in a trapped system. 

The motivation behind this work is that we believe the 
theory of 2D- YBF can be used in studying the ground 
state property of the strongly correlated liquid formed 
at the melted phase of many vortex system in the dilute 
trapped Bose gas [1, ■ In the latter system Bose- 

Einstein condensate (BEC) is completely depleted by 
quantum fluctuations, and quantum liquids appear with 
excitations that can carry fractional statistics. In this 
situation vortices can be treated as fundamental objects 
forming themselves as strongly correlated quantum fluid. 
This idea has been used in explaining certain features 
of Hall conductance and magnetization experiments in 
high-T c superconductors 0, 0, EH . On the other hand, 
the repulsive nature of vortices interacting in the short 
ranges is well described by the soft-core logarithmic po- 
tential. At the edge of the trap the potential is smoothed 
out by exp(—x)/y/i that decays to zero making it com- 
putational tractable and rule out any instability that 
may occur at the trap boundary. These characteristic 
making it plausible for the vortex liquid to be mapped 
as Yukawa bosons in a similar fashion to one adapted by 
Nelson and Seung Q for the YBF-FLL in the high T c 
superconductors. 

The paper is outlined as follows. In scction|TT]we intro- 
duce the Density functional Theory formalism in showing 
that by choosing a particular type of auxiliary external 
trapping potential the non-interacting system of bosons 
can be mapped to a real interacting system. We present 
the numerical calculation the Kohn-Sham equation and 
Gross-Pitacvskii (GP) in section [TTTT Then the numerical 
results will be discussed and concluded. 
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V ext + V int (1) 



where Hq = J drip'(r) 



A-V 2 - M J 1&(r) and V(\v - 
r'|) is the inter-atomic interaction potential. Here V ext (r) 
is the external trapping potential while m and \i are the 
atomic mass and the chemical potential respectively. The 
annihilation and creation field operators are denoted by 
■0 (r) and ip(r') respectively and obey Bose-Einstein com- 
mutation relations: 



[ip(r),^{r')] = 5(r-r') 
[V;(r),0(r')] = [0 t (r) 1 0t( r ')] =O . 



(2) 



Let us denote the ground state of the system as | ^> ) so 
the ground state energy is defined as E a = ($ | H \^ ) 
and the ground state density by n (r) = (\P | tp^tp \^ ). 
The Hohenberg-Kohn (HK) theorem [lj| guarantees that 
there exists a unique functional of the density, 



F[n}=H [n} + V mt [n], 



(3) 



irrespective of the external potential. The theorem 
was originally proved for fermions but its generalization 
also covers bosons. Following HK, we can write the total 
energy functional of the system as following, 



E[n] = F[n] + J drV ext (r)n(r) . 



(4) 



Determination of the ground state energy E Q follows 
by imposing the stationary conditions 



II. APPLICATION OF DFT THEORY IN THE 
YUKAWA BOSONS SYSTEM 



The Density functional theory (DFT) which is origi- 
nally based on the notion that for a many-electron sys- 
tem there is a one-to-one mapping between the external 
potential and the electron density: v ext (r) <-> p(r). In 
other words, the density is uniquely determined given a 
potential, and vice versa. All properties are therefore a 
functional of the density, because the density determines 
the potential, which determines the Hamiltonian, which 
determines the energy and the wave function. Following 
this train of reasoning, the inhomogeneous dilute system 
of N interacting bosons can be described within the sec- 
ond quantization language as 



SE[n] 
Sn(r) 







(5) 



where we will obtain the ground state density n (r) 
that is uniquely determined by the choice of our external 
potential V ex t- In general, the Hohenberg-Kohn theo- 
rem does not provide us with a computational scheme 
to determine the ground state energy. This is provided 
by the Kohn-Sham (KS) procedure [l5[. The idea is to 
use an auxiliary system (non-interacting reference sys- 
tem) and look for and external potential V* xt such that 
the noninteracting has the same ground state density as 
the real, interacting system. We write the Hamiltonian 
of the auxiliary system in the following form: 



H s = J dr^(r 

+ [ dv^{v)Vi xt {v)^{v). 



h2 n2 
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V>(r) 
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Thus for the ground state |W G ) we can define the unique 
total ground state energy functional of the auxiliary can 
be written as 



E s [n s ] = F s [n s ] + J dr n s (r)V e s xt (r) 



(7) 



By the KS scheme we note that E s [n s ] can be approx- 
imated to E[n] or in other word the density of the auxil- 
iary system n s (r) is equivalent to the real system n(r) by 
choosing a proper choice of auxiliary external potential 
V s t 

v ext 

Using the above argument and comparing the energy 
term F s [n s ] of Eq. (J7J) with Eqs. $5§, we can deduce the 
following functional relation: 



F[n} = F s [n}+V H [n}+F xc [n] 



(8) 



where the second term Vh is called the Hartree-energy 
defined as 



V, 



H 



l 1 



dr'V(|r-r'|)n(r)n(r') 



(9) 



The last term in Eq. ((SJ) represents the exchange- 
correlation energy E xc [n(r)] that includes all the con- 
tributions to the interaction energy beyond mean field 
Hatree term. Calculating the variational derivatives in 
Eq. © using Eq.©, one finds 



f^ + Mr) + ^l + Ur)=fl (10) 
dn(r) on(r) 

where the Hartree field read 

V H (v) = J dr'V(\v-v'\)n(r'). (11) 

Performing a similar variational calculation on Eq. ([7]) 
we deduce that the density of the auxiliary system is 
identical to the actual system if 



V^ t (v) = V ext (v) + V H (v) + V xc (v), 



(12) 



where we have introduced in Eq (fT2"|) the exchange- 
correlation potential V xc = SF xc [n(r)]/Sn(r) which is un- 
known for most of the system of interest and thus one 
has to resort to approximations such as the Local Den- 
sity Approximation (LDA) . This will be dealt with in the 
following subsequent section. 



A. The local density approximation 



such approximation to be suggested was the Local Den- 
sity Approximation, or LDA. The idea behind the LDA 
is very simple; it just ignores the non-local aspects of the 
functional dependence of V xc (r). The true form of V xc (r) 
will depend not only on the local density n(r) but also on 
n at all other points r' and this functional dependence is 
in general not known. This difficulty is avoided with the 
assumption that V xc depends only on the local density 
n(r), and that -E^n] can thus be written as 



E xc [n{r) 



dvE^ m [p]\ P - 



(13) 



where E x ° m [p] = pe xc [p] and e xc [p] is the exchange 
correlation energy of a homogeneous system with uni- 
form density p. The Functional derivatives of the above 
relation reads : 



V xc [n{l)} = 



5E xc [n(r)] d(pe xc [p\) 



Sn(r) 



dp 



\p— m(r) 



(14) 



The available numerical data of Magro and Ceperley 
[H and Strepparola et. al 0] permit one to obtain in- 
formation on the homogeneous excess free energy f ex [n] 
rather than the homogeneous exchange correlation en- 
ergy /xcM ■ Thus it is much convenient to work with 
excess free energy defined as [l6| 



F ex [n{r)} = V H [n(r)} + F xc [n(r)} ■ (15) 

In general, the excess correlation functional energy 
F ea; [n(r)] in a DFT calculation is not known exactly. One 
can resort to approximations such as the Local density 
approximation (LDA) which reads, 



J?e*[n(r)]« / dri^r>]U„ (r) 



(16) 



where F^° m [n] — nf^™ 1 [n] . The functional derivatives 
(excess-correlation potential) of the above relation can be 
written as: 



V ex (r,n(r)) = 



SF ex [n(r)} d(nf ex [n]) 



6n(r) 



dn 



i — m(r) • 



(17) 



Information of the homogeneous excess correlation en- 
ergy fex[n] can be obtained by subtracting the kinetic 
energy from the total ground state energy of a homo- 
geneous system with N boson. In Fig. ((2]) a plot of the 
DMC data of Magro and Ceperley 0] and the STLS data 
of Strepparola et. al. Q are depicted along with the fol- 
lowing fit: 



Before we can actually implement the Kohn-Sham for- 
malism, we have to devise some workable approximation 
for the exchange-correlation potential V xc (r). The first 



fex[p] = - 



_2p_ 

log(p) 



(18) 



4 



0.045 




5 



parola et. al. (crosses) compared to the fit function Eq. 
(fT8)> (dashed line). 



III. NUMERICAL CALCULATION FOR THE 
KOHN-SHAM EQUATION. 

We summarize here the algorithm in our numerical 
calculations which are commonly used in the literature 
to minimize energy functional E[fy, "J*] =< ^\H + 
V e s ;ct | v E' > where vff = \Jn(r) exp(i(p) in which the phase 
<f> fixes the velocity of the fluid through the relation 
v = \/<j)/m . We minimize .Ef 1 ]/, \&*] by assuming the 
normalization condition J dr^*^ = N, to obtain the 
non-linear the time independent non-linear Kohn-Sham 
equation : 



2n(x) 



log[rc(x)] 



*S>(x) = . 



(19) 



In the above equation we have incorporated the auxil- 
iary external potential based on Eq. (fT2")) and Eqs. (fT5|) - 
(|18[) using an isotropic planar trapping external potential 
V ex t = l/2mu! 2 r 2 with u) as radial frequency. We have 
scaled all length and energy by harmonic oscillator unit 
a/j = \J h/muj and huj/2 respectively. We obtained the 
ground state solution of the system by numerical iter- 
ation method by using a two-step Crank-Nicholson dis- 
cretization technique. 

The density profiles as a function of N is shown in Fig. 
([3]) . The profile shows a maximum at the center of trap 
and it decreases monotonically with x. To make a quan- 
titative measurement, we fix the solution for Eq. (|19|) for 
N = 100 (ni) as a standard and compare the differences 
in density for the other two values at the trap center 
(An = n 2 (0) - m(0) and An = 713(0) - ni(0) ). Here 
112(0) and 712(0) corresponds to the solution of Eq. (fTi?|) 
for N = 150 and N = 200 respectively at the center of 
trap. We found the differences An and An decreases with 
An/ni and An/711 by as large as 17% to 26%. The ex- 
planation of these results is quite straight forward. The 
areal integral of the density n(r) yields the total number 



FIG. 3: Top panel : Density profiles of the YBF n(x) (in 
arbitrary units) as a function of x (in ^Jh/muj units ) for 
various number of atoms N = 200 (solid line), 150 (dashed 
line) and 100 (dashed dot line). 



of atoms in the system. Hence larger number of atoms 
produced a much profound profile. The central density 
of the cloud decreases rapidly with increasing N but the 
density distribution is flattened due to stronger repulsion 
between particles. 

We also compare the YBF model with the conventional 
Gross-Piteavskii (GP) equation. This can be done by 
taking a contact potential gS(\r — r'\) in Eq. (fTTj) . where 
g = 47r/7 2 /mlog 1 2 H 2 /(m/ia 2 ) | is the coupling parameter 
for the two-dimensional BEC [13, [H, El incorporating 
the s-wave scattering legth a. We assumed a negligi- 
ble exchange correlation between the atoms (V xc = 0) 
and would like to stress here that the condensate den- 
sity is approximated to the total density ncp(r) at abso- 
lute zero temperature (neglecting quantum fluctuation). 
Based on the similar arguments in obtaining Eq. (|19[) , we 
obtain the time independent Gross-Pitaevskii equation : 



2 



+ x 2 + griGp(x) 



$?(x) = (j,$(x) , 



(20) 



where g is the dimensionless scaled coupling parame- 
ter. We have also calculated the total energy of YBF and 
GP model as a function of the number of atoms and the 
result is shown in Figure Increasing N we observe 
an increase of both interaction and harmonic oscillator 
potential energy for both GP and YBF models. The lat- 
ter effect follows from the expansion of the cloud. On 
the contrary, the kinetic energy per particle decreases 
because the density profile is flattened. The GP energy 
curve (dash-line) remains well above YBF curve (solid 
line) for all ranges of N. The energetic superiority of 
GP solution demonstrate the strongly repulsive nature 
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scheme in reproducing the result of Monte-Carlo simula- 
tion of Magro and Ceperley [3] for the liquid phase in 
a harmonic trap. Physically sensible result through the 
Local Density Approximation is obtained for the trapped 
system by knowing the homogeneous exchange correla- 
tion energy. The ground state properties (density pro- 
files and total energies) has been obtained by solving the 
non-linear equations Eqs. (TT9)) and . Our results 
shows that bosons interacting with Yukawa potential is 
energetically favourable compared to the contact delta 
potential for all ranges of N considered in this work. 
The results have so far not been verified. In light of 
the above, we trust that the YBF model through an ap- 
propriate mapping (boson- vortex) can be efficiently used 
to study strongly correlated liquid formed at the melted 
phase of vortices in a harmonically trapped rotating Bosc 
gases. 



FIG. 4: Total energies (in hui/2 units) as a function of N. 

of delta potential compared to the soft-core nature of 
Yukawa potential as the number of atoms increases in 
the system. 

In summary, we have studied the system of N Yukawa 
bosons in a 2D harmonic trap at absolute zero temper- 
ature. The central issue of this work is the use of Den- 
sity Functional Theory formalism within the Kohn-Sham 
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